# Medicaid expansion and political participation
# J.R. Jolin, N.M. Kavanagh
# December 10, 2022

# Please direct questions about this script file to nolankavanagh@fas.harvard.edu.

# Clear R environment
rm(list = ls())

# Load packages
library(dplyr)        # Analysis tools
library(tidyverse)    # Analysis tools
library(naniar)       # Analysis tools
library(stringr)      # Analysis tools
library(readstata13)  # Analysis tools
library(fastDummies)  # Analysis tools
library(lfe)          # Modeling tools
library(survey)       # Modeling tools
library(estimatr)     # Modeling tools
library(ggplot2)      # Graphing tools
library(gridExtra)    # Graphing tools
library(cowplot)      # Graphing tools
library(usmap)        # Mapping tools
library(modelsummary) # Table tools
library(kableExtra)   # Table tools
library(arsenal)      # Table tools
library(utils)        # Table tools

# Set working directory
#setwd("/Users/jamesjolin/Desktop/Dropbox/Medicaid and political participation/")
#setwd("/Users/nolankavanagh/Dropbox/Medicaid and political participation/")

##############################################################################
# Preparing CCES data set
##############################################################################

# Read in data
cces <- readRDS("CES Common Content 2006-2021/cumulative_2006-2021.rds")

# Subset to citizens only
# cces <- subset(cces, cces$citizen == 1)

# Age
cces <- cces %>% mutate(
  age_r = case_when(
    age %in% c(18:30)  ~ "18 to 30",
    age %in% c(31:40)  ~ "31 to 40",
    age %in% c(41:50)  ~ "41 to 50",
    age %in% c(51:60)  ~ "51 to 60",
    age %in% c(61:70)  ~ "61 to 70",
    age %in% c(71:150) ~ "71 and older"))

# Gender
cces <- cces %>% mutate(
  gender_r = case_when(
    gender %in% c(1) ~ "Male",
    gender %in% c(2) ~ "Female"))
cces$gender_r <- factor(cces$gender_r, levels = c("Male", "Female"))

# Race/ethnicity
cces <- cces %>% mutate(
  race_r = case_when(
    race %in% c(1)   ~ "White, non-Hispanic/Latino",
    race %in% c(2)   ~ "Black, non-Hispanic/Latino",
    race %in% c(3)   ~ "Hispanic/Latino",
    race %in% c(4)   ~ "Asian",
    race %in% c(5)   ~ "Native American",
    race %in% c(6:8) ~ "Other"))
cces$race_r <- factor(cces$race_r, levels = c("White, non-Hispanic/Latino", "Black, non-Hispanic/Latino", "Hispanic/Latino", "Asian", "Native American", "Other"))

# White/non-white
cces <- cces %>% mutate(
  nonwhite_r = case_when(
    race == 1   ~ 0,
    race %in% c(2:8) ~ 1))

# Education
cces <- cces %>% mutate(
  educ_r = case_when(
    educ %in% c(1)   ~ "Less than high school",
    educ %in% c(2)   ~ "High school graduate",
    educ %in% c(3:4) ~ "Some college/2-year degree",
    educ %in% c(5)   ~ "College graduate",
    educ %in% c(6)   ~ "Postgraduate degree"))
cces$educ_r <- factor(cces$educ_r, levels = c("Less than high school", "High school graduate", "Some college/2-year degree", "College graduate", "Postgraduate degree"))

# Marital status
cces <- cces %>% mutate(
  marstat_r = case_when(
    marstat %in% c(1,6)   ~ "Married or Domestic Partnership",
    marstat == 4          ~ "Widowed",
    marstat %in% c(2,3,5) ~ "Single / Separated / Divorced"))
cces$`marstat_r` <- factor(cces$`marstat_r`, levels = c("Married or Domestic Partnership", "Widowed", "Single / Separated / Divorced"))

# Has a child
cces <- cces %>% mutate(
  has_child_r = case_when(
    has_child == "Yes" ~ 1,
    has_child == "No"  ~ 0))

# Employment status
cces <- cces %>% mutate(
  employ_r = case_when(
    employ %in% c("Full-Time", "Part-Time")             ~ "Employed",
    employ %in% c("Temporarily Laid Off", "Unemployed") ~ "Unemployed",
    employ %in% c("Retired")                            ~ "Retired",
    employ %in% c("Permanently Disabled")               ~ "Disabled",
    employ %in% c("Homemaker", "Student", "Other")      ~ "Other"
  ))
cces$employ_r <- factor(cces$employ_r, levels = c("Employed", "Unemployed", "Retired", "Disabled", "Other"))

# Family income
cces <- cces %>% mutate(
  income_r = case_when(
    faminc %in% c("Less than 10k", "10k - 20k")   ~ "$19,999 or less",
    faminc %in% c("20k - 30k", "30k - 40k")       ~ "$20,000 to $39,999",
    faminc %in% c("40k - 50k", "50k - 60k")       ~ "$40,000 to $59,999",
    faminc %in% c("60k - 70k", "70k - 80k")       ~ "$60,000 to $79,999",
    faminc %in% c("80k - 100k")                   ~ "$80,000 to $99,999",
    faminc %in% c("100k - 120k", "120k - 150k", "150k+")  ~ "$100,000 or more",
    faminc %in% c("Prefer not to say", "Skipped") ~ "Prefer not to say"
  ))
cces$income_r <- factor(cces$income_r, levels = c("$19,999 or less", "$20,000 to $39,999", "$40,000 to $59,999", "$60,000 to $79,999", "$80,000 to $99,999", "$100,000 or more", "Prefer not to say"))

# Union status
cces <- cces %>% mutate(
  union_r = case_when(
    union == 1 ~ "Yes, currently",
    union == 2 ~ "Yes, formerly",
    union == 3 ~ "No, never"
  ))

# Low-income households
cces <- cces %>% mutate(
  low_income = case_when(
    faminc %in% c("Less than 10k", "10k - 20k", "20k - 30k") ~ 1,
    faminc %in% c("30k - 40k", "40k - 50k", "50k - 60k", "60k - 70k", "70k - 80k",
                  "80k - 100k", "100k - 120k", "120k - 150k", "150k+") ~ 0))

# Define expansion status
cces <- cces %>% mutate(
  expansion = case_when(
    st %in% c("AR","AZ","CA","CO","CT","DC","DE","HI","IA","IL","KY","MA","MD",
              "MI","MN","NM","ND","NH","NJ","NY","OH","NV","OR","RI","VT","WA","WV") ~ 1,
    st %in% c("AL","FL","GA","KS","MS","NC","SC","TN","TX","WI","WY") ~ 0,
    st %in% c("ID","MO","NE","OK","SD","UT") ~ 0)) # Expanded in 2020 or later
# Exclude: AK (2015), IN (2015), LA (2016), MT (2016), PA (2015), VA (2019), ME (2019)

# Define expansion label
cces <- cces %>% mutate(
  `Expansion status` = case_when(
    expansion == 0 ~ "Non-expansion",
    expansion == 1 ~ "Expansion"))
cces$`Expansion status` <- factor(cces$`Expansion status`, levels = c("Expansion", "Non-expansion"))

# Define expansion years
cces <- cces %>% mutate(
  time = case_when(
    year %in% c(2010:2013) ~ 0,
    year %in% c(2014:2019) ~ 1))

# Identify as Democrats
cces <- cces %>% mutate(
  democrat = case_when(
    pid3 %in% c(1)   ~ 1,
    pid3 %in% c(2:5) ~ 0))

# Identify as Republicans
cces <- cces %>% mutate(
  republican = case_when(
    pid3 %in% c(2)     ~ 1,
    pid3 %in% c(1,3:5) ~ 0))

# Ideology self-placement
cces <- cces %>% mutate(
  ideo5_r = case_when(
    ideo5 == "Very Liberal"      ~  2,
    ideo5 == "Liberal"           ~  1,
    ideo5 == "Moderate"          ~  0,
    ideo5 == "Conservative"      ~ -1,
    ideo5 == "Very Conservative" ~ -2))

# Partisan strength
cces <- cces %>% mutate(
  pid7_r = case_when(
    pid7 == 1 ~  3,
    pid7 == 2 ~  2,
    pid7 == 3 ~  1,
    pid7 == 4 ~  0,
    pid7 == 5 ~ -1,
    pid7 == 6 ~ -2,
    pid7 == 7 ~ -3))

# Presidential approval
cces <- cces %>% mutate(
  approval_pres_r = case_when(
    approval_pres == 1 ~ 1,
    approval_pres == 2 ~ 0.5,
    approval_pres %in% c(5,6) ~ 0,
    approval_pres == 3 ~ -0.5,
    approval_pres == 4 ~ -1))

# Representative approval
cces <- cces %>% mutate(
  approval_rep_r = case_when(
    approval_rep == "Strongly Approve" ~ 1,
    approval_rep == "Approve / Somewhat Approve" ~ 0.5,
    approval_rep %in% c("Never Heard / Not Sure","Never Heard of Person","Never Heard of this Person","Neither Approve nor Disapprove") ~ 0,
    approval_rep == "Disapprove / Somewhat Disapprove" ~ -0.5,
    approval_rep == "Strongly Disapprove" ~ -1))

# Senator 1 approval
cces <- cces %>% mutate(
  approval_sen1_r = case_when(
    approval_sen1 == "Strongly Approve" ~ 1,
    approval_sen1 == "Approve / Somewhat Approve" ~ 0.5,
    approval_sen1 %in% c("Never Heard / Not Sure","Never Heard of Person","Never Heard of this Person","Neither Approve nor Disapprove") ~ 0,
    approval_sen1 == "Disapprove / Somewhat Disapprove" ~ -0.5,
    approval_sen1 == "Strongly Disapprove" ~ -1))

# Senator 2 approval
cces <- cces %>% mutate(
  approval_sen2_r = case_when(
    approval_sen2 == "Strongly Approve" ~ 1,
    approval_sen2 == "Approve / Somewhat Approve" ~ 0.5,
    approval_sen2 %in% c("Never Heard / Not Sure","Never Heard of Person","Never Heard of this Person","Neither Approve nor Disapprove") ~ 0,
    approval_sen2 == "Disapprove / Somewhat Disapprove" ~ -0.5,
    approval_sen2 == "Strongly Disapprove" ~ -1))

# Average approval for senators
cces <- cces %>% mutate(
  approval_sen_r = (approval_sen1_r + approval_sen2_r)/2)

# Subset to 2010-2019
cces <- cces %>% filter(
  year %in% c(2010:2019))

# Subset to complete cases for covariates
cces <- cces %>% filter_at(
  vars(age_r, gender_r, race_r, educ_r, employ_r, income_r,
       marstat_r, has_child_r, union_r), all_vars(!is.na(.)))

##############################################################################
# Demographics of respondents
##############################################################################

# Create summary statistics -- unweighted
demo_table <- summary(tableby(~ gender_r + race_r + income_r + race_r + employ_r + marstat_r + has_child_r + union_r, cces))

# Created summary statistics -- weighted
demo_table_w <- summary(tableby(~ gender_r + race_r + income_r + race_r + employ_r + marstat_r + has_child_r + union_r, cces, weights = cces$weight_cumulative))

# Export summary tables
demotables <- knitr::kable(demo_table, format="latex")
writeLines(demotables,"demographicinfo.tex")
demotables <- knitr::kable(demo_table_w, format="latex")
writeLines(demotables,"demographicinfoweighted.tex")

##############################################################################
# Map of included states
##############################################################################

# Generate dataframe of states
states_df <- cces %>%
  group_by(st) %>%
  summarise(expansion = median(expansion))
states_df <- as.data.frame(states_df)
colnames(states_df) <- c("state", "expansion")

# Treat expansion as factor
states_df$expansion <- factor(states_df$expansion, levels=c(1,0))

# Map for included states
map_medicaid <- plot_usmap(regions="states", data=states_df,
                           values="expansion", size=0.4) +
  theme(legend.position = "right",
        text = element_text(size=10, face="bold")) +
  scale_fill_manual(name   = "Expansion status",
                    values = c("#00BFC4","#F8766D", "dark grey"),
                    labels = c("Expanded in 2014",
                               "Didn't expand during 2014-2019",
                               "Not included (expanded during 2015-2019)"))

# Export figure
ggsave(plot=map_medicaid, file="Map of states.pdf", width=8, height=4, units='in', dpi=600)

##############################################################################
# Custom functions for subsequent models and graphics
##############################################################################

# Custom function to create all subsequent regression tables, for efficiency of code; to present other coefficients, use modelsummary()
createresults <- function(models_to_tabulate, title_enter, name_file) {
  return(modelsummary(models_to_tabulate,
                      gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
                      coef_omit   = "^(?!expansion:time)",
                      title       = title_enter,
                      coef_rename = c("expansion:time" = "Expansion",
                                      "expansion:time:low_income" =
                                        "Expansion × Low-Income",
                                      "expansion:time:nonwhite_r" =
                                        "Expansion × Non-White"),
                      add_rows    = cov_row,
                      stars       = c("*"=0.05, "**"=0.01, "***"=0.001),
                      output      = name_file))
}

##############################################################################
# Difference-in-difference models: Party identification
##############################################################################

# Models for percent Democrat
dem_1 <- felm(democrat ~ expansion*time
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
dem_2 <- felm(democrat ~ expansion*time +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
dem_3 <- felm(democrat ~ expansion*time*low_income +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
dem_4 <- felm(democrat ~ expansion*time*nonwhite_r +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)

# Models for percent Republican
rep_1 <- felm(republican ~ expansion*time
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
rep_2 <- felm(republican ~ expansion*time +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
rep_3 <- felm(republican ~ expansion*time*low_income +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
rep_4 <- felm(republican ~ expansion*time*nonwhite_r +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)

# Additional row in table
cov_row <- as.data.frame(
  rbind(cbind("Covariates",          t(rep(c("No",  "Yes", "Yes", "Yes"), 1))),
        cbind("State fixed effects", t(rep("Yes", 4))),
        cbind("Year fixed effects",  t(rep("Yes", 4))),
        cbind("Cluster robust SEs",  t(rep("Yes", 4)))))

# Present regression results
createresults(list("(1)" = dem_1,
                   "(2)" = dem_2,
                   "(3)" = dem_3,
                   "(4)" = dem_4),
              "Proportion of respondents identifying as Democrats after Medicaid expansion", "dem_id_table.tex")

createresults(list("(1)" = rep_1,
                   "(2)" = rep_2,
                   "(3)" = rep_3,
                   "(4)" = rep_4),
              "Proportion of respondents identifying as Republicans after Medicaid expansion",
              "rep_id_table.tex")

##############################################################################
# Graph: Party identification
##############################################################################

# Define study designs
design <- svydesign(id=~state, weights=~weight_cumulative, data=cces)

# Tables for percent parties
table <- svyby(~democrat + republican, ~year + `Expansion status`,
               design, svymean, na.rm=TRUE, vartype="ci")

# Treat years as numeric
table$year <- as.numeric(table$year)

# Unadjusted plots: Percent Democrat
plot_dem <- ggplot(table, aes(x=year, y=democrat,
                              ymin  = ci_l.democrat,
                              ymax  = ci_u.democrat,
                              group =`Expansion status`,
                              fill  =`Expansion status`,
                              color =`Expansion status`)) +
  geom_vline(xintercept=2013.75, color="black", linetype="dashed") +
  geom_ribbon(alpha=0.25, linetype = 0) +
  geom_point() + geom_line() +
  theme_test() +
  theme(legend.position    = "none",
        text               = element_text(size = 10, face = "bold"),
        axis.ticks         = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25),
        panel.grid.major.y = element_line(color = "gray", size = 0.25),
        panel.grid.minor.y = element_line(color = "light gray", size = 0.25)) +
  xlab("Year") + ylab("Percent identifying\nas Democrats") +
  annotate("text", x=2014.25, y=0.22, label="Medicaid expansion",
           fontface=2, size=3, hjust=0, vjust=0) +
  scale_x_continuous(limits = c(2009.5, 2019.5),
                     breaks = seq(2010, 2019, 2)) +
  scale_y_continuous(limits = c(0.15,0.45),
                     labels = scales::percent_format(accuracy=1),
                     breaks = c(0.2,0.3,0.4),
                     minor_breaks=c(0.15,0.25,0.35,0.45)) +
  scale_color_manual(values = c("#00BFC4","#F8766D")) +
  scale_fill_manual(values = c("#00BFC4","#F8766D"))

# Unadjusted plots: Percent Republican
plot_rep <- ggplot(table, aes(x=year, y=republican,
                              ymin  = ci_l.republican,
                              ymax  = ci_u.republican,
                              group =`Expansion status`,
                              fill  =`Expansion status`,
                              color =`Expansion status`)) +
  geom_vline(xintercept=2013.75, color="black", linetype="dashed") +
  geom_ribbon(alpha=0.25, linetype = 0) +
  geom_point() + geom_line() +
  theme_test() +
  theme(legend.position    = "none",
        text               = element_text(size = 10, face = "bold"),
        axis.ticks         = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25),
        panel.grid.major.y = element_line(color = "gray", size = 0.25),
        panel.grid.minor.y = element_line(color = "light gray", size = 0.25)) +
  xlab("Year") + ylab("Percent identifying\nas Republicans") +
  annotate("text", x=2014.25, y=0.38, label="Medicaid expansion",
           fontface=2, size=3, hjust=0, vjust=0.5) +
  scale_x_continuous(limits = c(2009.5, 2019.5),
                     breaks = seq(2010, 2019, 2)) +
  scale_y_continuous(limits = c(0.15,0.45),
                     labels = scales::percent_format(accuracy=1),
                     breaks = c(0.2,0.3,0.4),
                     minor_breaks=c(0.15,0.25,0.35,0.45)) +
  scale_color_manual(values = c("#00BFC4","#F8766D")) +
  scale_fill_manual(values = c("#00BFC4","#F8766D"))

# Legend
legend <- get_legend(
  ggplot(table, aes(x=year, y=democrat,
                    ymin  = ci_l.democrat,
                    ymax  = ci_u.democrat,
                    group =`Expansion status`,
                    fill  =`Expansion status`,
                    color =`Expansion status`)) +
    geom_ribbon(alpha=0.25, linetype = 0) +
    geom_point() + geom_line() +
    theme(legend.position = "right",
          text = element_text(size = 10, face = "bold")) +
    scale_color_manual(values = c("#00BFC4","#F8766D")) +
    scale_fill_manual(values = c("#00BFC4","#F8766D")))

##############################################################################
# Difference-in-difference models: Ideology self-placement
##############################################################################

# Models for ideology
ideo_1 <- felm(ideo5_r ~ expansion*time
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
ideo_2 <- felm(ideo5_r ~ expansion*time +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
ideo_3 <- felm(ideo5_r ~ expansion*time*low_income +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
ideo_4 <- felm(ideo5_r ~ expansion*time*nonwhite_r +
                 age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
               | st + year | 0 | st, cces, weights=cces$weight_cumulative)

# Additional row in table
cov_row <- as.data.frame(
  rbind(cbind("Covariates",          t(c("No", "Yes", "Yes", "Yes"))),
        cbind("State fixed effects", t(rep("Yes", 4))),
        cbind("Year fixed effects",  t(rep("Yes", 4))),
        cbind("Cluster robust SEs",  t(rep("Yes", 4)))))

# Present regression results
createresults(list("(1)" = ideo_1,
                   "(2)" = ideo_2,
                   "(3)" = ideo_3,
                   "(4)" = ideo_4),
              "Mean ideology self-placement after Medicaid expansion",
              "ideology_table.tex")

##############################################################################
# Graph: Ideology self-placement
##############################################################################

# Table for ideology
table <- svyby(~ideo5_r, ~year + `Expansion status`,
               design, svymean, na.rm=TRUE, vartype="ci")

# Treat years as numeric
table$year <- as.numeric(table$year)

# Unadjusted plots
plot_ideo <- ggplot(table, aes(x=year, y=ideo5_r,
                              ymin  = ci_l,
                              ymax  = ci_u,
                              group =`Expansion status`,
                              fill  =`Expansion status`,
                              color =`Expansion status`)) +
  geom_vline(xintercept=2013.75, color="black", linetype="dashed") +
  geom_ribbon(alpha=0.25, linetype = 0) +
  geom_point() + geom_line() +
  theme_test() +
  theme(legend.position    = "none",
        text               = element_text(size = 10, face = "bold"),
        axis.ticks         = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25),
        panel.grid.major.y = element_line(color = "gray", size = 0.25),
        panel.grid.minor.y = element_line(color = "light gray", size = 0.25)) +
  xlab("Year") + ylab("Mean ideology self-placement\n(-2 cons. to 2 lib.)") +
  annotate("text", x=2014.25, y=0.187, label="Medicaid expansion",
           fontface=2, size=3, hjust=0, vjust=0.5) +
  scale_x_continuous(limits = c(2009.5, 2019.5),
                     breaks = seq(2010, 2019, 2)) +
  scale_y_continuous(limits = c(-0.55,0.3),
                     breaks = seq(-0.5, 0.25, by=0.25)) +
  scale_color_manual(values = c("#00BFC4","#F8766D")) +
  scale_fill_manual(values = c("#00BFC4","#F8766D"))

##############################################################################
# Difference-in-difference models: Partisan strength
##############################################################################

# Models for partisan strength
pids_1 <- felm(pid7_r ~ expansion*time
               | st + year | 0 | st, cces, weights=cces$weight_cumulative)
pids_2 <- felm(pid7_r ~ expansion*time +
                 age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
               | st + year | 0 | st, cces, weights=cces$weight_cumulative)
pids_3 <- felm(pid7_r ~ expansion*time*low_income +
                 age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
               | st + year | 0 | st, cces, weights=cces$weight_cumulative)
pids_4 <- felm(pid7_r ~ expansion*time*nonwhite_r +
                 age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
               | st + year | 0 | st, cces, weights=cces$weight_cumulative)

# Additional row in table
cov_row <- as.data.frame(
  rbind(cbind("Covariates",          t(c("No",  "Yes", "Yes", "Yes"))),
        cbind("State fixed effects", t(rep("Yes", 4))),
        cbind("Year fixed effects",  t(rep("Yes", 4))),
        cbind("Cluster robust SEs",  t(rep("Yes", 4)))))

# Present regression results
createresults(list("(1)" = pids_1,
                   "(2)" = pids_2,
                   "(3)" = pids_3,
                   "(4)" = pids_4),
              "Mean partisanship strength after Medicaid expansion",
              "pid_strength_table.tex")

##############################################################################
# Graph: Partisan strength
##############################################################################

# Tables for partisan strength
table <- svyby(~pid7_r, ~year + `Expansion status`,
               design, svymean, na.rm=TRUE, vartype="ci")

# Treat years as numeric
table$year <- as.numeric(table$year)

# Unadjusted plot
plot_pid7 <- ggplot(table, aes(x=year, y=pid7_r,
                               ymin  = ci_l,
                               ymax  = ci_u,
                               group =`Expansion status`,
                               fill  =`Expansion status`,
                               color =`Expansion status`)) +
  geom_vline(xintercept=2013.75, color="black", linetype="dashed") +
  geom_ribbon(alpha=0.25, linetype = 0) +
  geom_point() + geom_line() +
  theme_test() +
  theme(legend.position    = "none",
        text               = element_text(size = 10, face = "bold"),
        axis.ticks         = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25),
        panel.grid.major.y = element_line(color = "gray", size = 0.25),
        panel.grid.minor.y = element_line(color = "light gray", size = 0.25)) +
  xlab("Year") + ylab("Mean partisan strength\n(-3 Rep. to 3 Dem.)") +
  annotate("text", x=2014.25, y=-0.3, label="Medicaid expansion",
           fontface=2, size=3, hjust=0, vjust=0.5) +
  scale_x_continuous(limits = c(2009.5, 2019.5),
                     breaks = seq(2010, 2019, 2)) +
  scale_y_continuous(limits = c(-0.45,0.85),
                     breaks = seq(-0.4, 0.8, by=0.4)) +
  scale_color_manual(values = c("#00BFC4","#F8766D")) +
  scale_fill_manual(values = c("#00BFC4","#F8766D"))

##############################################################################
# Graph: Combined partisanship and ideology panels
##############################################################################

# Combine all four panels
did_plots_combined <- plot_grid(plot_dem, plot_rep, legend, plot_pid7, plot_ideo, legend,
                       rel_widths = c(1, 1, 0.4), nrow = 2, vjust = 1,
                       labels=c("A.", "B.", "", "C.", "D.", ""))

# Export figure
ggsave(plot=did_plots_combined, file="Combined PID,ideo plots.pdf",
       width=8, height=5.5, units='in', dpi=600)

##############################################################################
# Difference-in-difference models: Approval of political leaders
##############################################################################

# Subset for Obama years
cces_pres_subset <- cces %>% filter(year %in% 2010:2016)

# Models for presidential approval
a_pres_1 <- felm(approval_pres_r ~ expansion*time
              | st + year | 0 | st, cces_pres_subset, weights=cces_pres_subset$weight_cumulative)
a_pres_2 <- felm(approval_pres_r ~ expansion*time +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces_pres_subset, weights=cces_pres_subset$weight_cumulative)
a_pres_3 <- felm(approval_pres_r ~ expansion*time*low_income +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces_pres_subset, weights=cces_pres_subset$weight_cumulative)
a_pres_4 <- felm(approval_pres_r ~ expansion*time*nonwhite_r +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces_pres_subset, weights=cces_pres_subset$weight_cumulative)

# Models for representative approval
a_rep_1 <- felm(approval_rep_r ~ expansion*time
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
a_rep_2 <- felm(approval_rep_r ~ expansion*time +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r 
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
a_rep_3 <- felm(approval_rep_r ~ expansion*time*low_income +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)
a_rep_4 <- felm(approval_rep_r ~ expansion*time*nonwhite_r +
                age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
              | st + year | 0 | st, cces, weights=cces$weight_cumulative)

# Models for mean senator approval
a_sen_1 <- felm(approval_sen_r ~ expansion*time
                | st + year | 0 | st, cces, weights=cces$weight_cumulative)
a_sen_2 <- felm(approval_sen_r ~ expansion*time +
                  age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r 
                | st + year | 0 | st, cces, weights=cces$weight_cumulative)
a_sen_3 <- felm(approval_sen_r ~ expansion*time*low_income +
                  age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
                | st + year | 0 | st, cces, weights=cces$weight_cumulative)
a_sen_4 <- felm(approval_sen_r ~ expansion*time*nonwhite_r +
                  age_r + gender_r + race_r + educ_r + employ_r + income_r + marstat_r + has_child_r + union_r
                | st + year | 0 | st, cces, weights=cces$weight_cumulative)

# Additional rows in table
cov_row <- as.data.frame(
  rbind(cbind("Covariates",          t(rep(c("No",  "Yes", "Yes", "Yes"),1))),
        cbind("State fixed effects", t(rep("Yes", 4))),
        cbind("Year fixed effects",  t(rep("Yes", 4))),
        cbind("Cluster robust SEs",  t(rep("Yes", 4)))))

# Present regression results
createresults(list("(1)" = a_pres_1, "(2)" = a_pres_2, "(3)" = a_pres_3, "(4)" = a_pres_4),
              "Mean approval of President Obama after Medicaid expansion",
              "pres_approval_table.tex")
createresults(list("(1)" = a_rep_1, "(2)" = a_rep_2, "(3)" = a_rep_3, "(4)" = a_rep_4),
              "Mean approval of respondent's representative after Medicaid expansion",
              "rep_approval_table.tex")
createresults(list("(1)" = a_sen_1, "(2)" = a_sen_2, "(3)" = a_sen_3, "(4)" = a_sen_4),
              "Mean approval of respondent's senators after Medicaid expansion","sen_approval_table.tex")

##############################################################################
# Graphs: Approval of officeholders
##############################################################################

# Table for approval
table <- svyby(~approval_pres_r + approval_rep_r + approval_sen_r,
               ~year + `Expansion status`,
               design, svymean, na.rm=TRUE, vartype="ci")

# Treat years as numeric
table$year <- as.numeric(table$year)

# Unadjusted plots
plot_pres <- ggplot(table, aes(x=year, y=approval_pres_r,
                              ymin  = ci_l.approval_pres_r,
                              ymax  = ci_u.approval_pres_r,
                              group =`Expansion status`,
                              fill  =`Expansion status`,
                              color =`Expansion status`)) +
  geom_vline(xintercept=2013.75, color="black", linetype="dashed") +
  geom_ribbon(alpha=0.25, linetype = 0) +
  geom_point() + geom_line() +
  theme_test() +
  theme(legend.position    = "none",
        text               = element_text(size = 10, face = "bold"),
        axis.ticks         = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25),
        panel.grid.major.y = element_line(color = "gray", size = 0.25),
        panel.grid.minor.y = element_line(color = "light gray", size = 0.25)) +
  xlab("Year") + ylab("Mean approval of Obama\n(-1 disapp. to 1 app.)") +
  annotate("text", x=2014.25, y=0.115, label="Medicaid expansion",
           fontface=2, size=3, hjust=0, vjust=0.5) +
  scale_x_continuous(breaks = seq(2010, 2019, 2), limits = c(2009.5, 2016.5)) +
  scale_y_continuous(breaks = c(-0.3, -0.15, 0, 0.15)) +
  coord_cartesian(xlim = c(2009.5, 2019.5), ylim = c(-0.3, 0.15)) +
  scale_color_manual(values = c("#00BFC4","#F8766D")) +
  scale_fill_manual(values = c("#00BFC4","#F8766D"))

plot_rep <- ggplot(table, aes(x=year, y=approval_rep_r,
                              ymin  = ci_l.approval_rep_r,
                              ymax  = ci_u.approval_rep_r,
                              group =`Expansion status`,
                              fill  =`Expansion status`,
                              color =`Expansion status`)) +
  geom_vline(xintercept=2013.75, color="black", linetype="dashed") +
  geom_ribbon(alpha=0.25, linetype = 0) +
  geom_point() + geom_line() +
  theme_test() +
  theme(legend.position    = "none",
        text               = element_text(size = 10, face = "bold"),
        axis.ticks         = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25),
        panel.grid.major.y = element_line(color = "gray", size = 0.25),
        panel.grid.minor.y = element_line(color = "light gray", size = 0.25)) +
  xlab("Year") + ylab("Mean approval of rep.\n(-1 disapp. to 1 app.)") +
  annotate("text", x=2014.25, y=-0.12, label="Medicaid expansion",
           fontface=2, size=3, hjust=0, vjust=0) +
  scale_x_continuous(limits = c(2009.5, 2019.5),
                     breaks = seq(2010, 2019, 2)) +
  scale_y_continuous(limits = c(-0.3, 0.15),
                     breaks = c(-0.3, -0.15, 0, 0.15)) +
  scale_color_manual(values = c("#00BFC4","#F8766D")) +
  scale_fill_manual(values = c("#00BFC4","#F8766D"))

plot_sen <- ggplot(table, aes(x=year, y=approval_sen_r,
                              ymin  = ci_l.approval_sen_r,
                              ymax  = ci_u.approval_sen_r,
                              group =`Expansion status`,
                              fill  =`Expansion status`,
                              color =`Expansion status`)) +
  geom_vline(xintercept=2013.75, color="black", linetype="dashed") +
  geom_ribbon(alpha=0.25, linetype = 0) +
  geom_point() + geom_line() +
  theme_test() +
  theme(legend.position    = "none",
        text               = element_text(size = 10, face = "bold"),
        axis.ticks         = element_blank(),
        panel.grid.major.x = element_line(color = "gray", size = 0.25),
        panel.grid.major.y = element_line(color = "gray", size = 0.25),
        panel.grid.minor.y = element_line(color = "light gray", size = 0.25)) +
  xlab("Year") + ylab("Mean approval of senators\n(-1 disapp. to 1 app.)") +
  annotate("text", x=2014.25, y=-0.185, label="Medicaid expansion",
           fontface=2, size=3, hjust=0, vjust=0.5) +
  scale_x_continuous(limits = c(2009.5, 2019.5),
                     breaks = seq(2010, 2019, 2)) +
  scale_y_continuous(limits = c(-0.3, 0.15),
                     breaks = c(-0.3, -0.15, 0, 0.15)) +
  scale_color_manual(values = c("#00BFC4","#F8766D")) +
  scale_fill_manual(values = c("#00BFC4","#F8766D"))

# Arrange plots
did_plots <- plot_grid(plot_pres, plot_rep, plot_sen, legend,
                       nrow = 2, vjust = 1, labels=c("A.", "B.", "C.", ""))

# Export figure
ggsave(plot=did_plots, file="Plots for approval.pdf",
       width=8, height=6, units='in', dpi=600)
